Large Scale Azimuthal Structures Of T\irbulence In Accretion Disks 

Dynamo triggered variability of accretion. 

M. Flock\ N. Dzyurkevich\ H. Klahr^ N. Turner^ ^Th. Henning' 
^Max Planck Institute for Astronomy, Konigstuhl 17, 691 17 Heidelberg, Germany 
^Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91 109, USA 



Received 



accepted 



-2- 
ABSTRACT 

We investigate the significance of large scale azimuthal, magnetic and velocity 
modes for the MRI turbulence in accretion disks. We perform 3D global ideal MHD 
simulations of global stratified proto-planetary disk models. Our domains span az- 
imuthal angles of 7r/4, 7r/2, n and In. We observe up to 100% stronger magnetic 
fields and stronger turbulence for the restricted azimuthal domain models n/l and nIA 
compared to the full In model. We show that for those models, the Maxwell Stress 
is larger due to strong axisymmetric magnetic fields, generated by the aO. dynamo. 
Large radial extended axisymmetric toroidal fields trigger temporal magnification of 
accretion stress. All models display a positive dynamo-a in the northern hemisphere 
(upper disk). The parity is distinct in each model and changes on timescales of 40 
local orbits. In model In, the toroidal field is mostly antisymmetric in respect to the 
midplane. The eddies of the MRI turbulence are highly anisotropic. The major wave- 
lengths of the turbulent velocity and magnetic fields are between one and two disk 
scale heights. At the midplane, we find magnetic tilt angles around 8-9° increasing 
up to 12 - 13° in the corona. We conclude that an azimuthal extent of n is sufficient to 
reproduce most turbulent properties in 3D global stratified simulations of magnetised 
accretion disks. 

Subject headings: accretion discs, magneto hydrodynamics (MHD), MHD Dynamo 
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1. Introduction 



Magneto-rotational instability (MRI) can generate MHD turbulenc e with an outward direct ed 



angular momentum transport driving accretion o nto the central object (iBalbus & Hawley 



Hawlev & Balbus 



1991 



Balbus & Hawlev 



19911: 



19981). A necessary condition is a good coupling 



between the gas and magnetic fields, e.g. a well-ionized gas. In proto-planetary disks, dust 



particles and low temperatures wil 



(Sanoetal 



20101 : 



2000 



Turner et al 



leming & Stone 



reduce the ionisation 



2003 



Inutsuka & Sano 



evel and therefor the MRI activity 



2005 



Wardle 2007 



Dzvurkevich et al 



20101) . Nevertheless, there are well-ionized regions with possible MRI activity. 



like the coronal region or the inn er or outer disk. Th e inner disk will be thermally ionized for 



temperatures greater then 1000 A' (lUmebayashi 



1983h. The outer disk wi 



1 be i onized by Cosmic 



20091). In our work we 



Rays for surface density values below 96g/cm^ (|Umebayashi & Nakand l! 
concentrate on well ionized disk regions. To model the evolution of proto-planetary disks and 
especially to describe the process of planet formation, we need to know detailed informations 
about the strength of the turbulence. Several processes, like the MHD dynamo or the toroidal 
field MRI, influence the turbulence level. The evolution of the magnetic and velocities fields at 
different scales has to be investigated. 



In the last decades, a large amou nt of local-box siniulations have been performed to 



study the small scale MRI turbulence (Brandenburg et al 



Matsumoto & Tajima 



seed magnetic fields ([Balbus & Hawley 



1995 



Stone et al. 



was analyzed through 



erauem & Paoaloizou 



(Gellert et al. 



2007 



inear calculations (jHawlev & Balbus 



1996 



1995 



Hawlev et al, 



1995, 



1996 



1996). The MRI works for both, vertical or toroidal 



1991 )■ The MRI launched with initial toroidal field 



Papaloizou & Terquem 



Riidiger et al. 



1992 



Foglizzo & Tagger 



19951: 



19971) and in Taylor-Couette experiments 



20071) . This experiments showed that most of the energy will 



be transported to th e m = 1 mode. A sim ilar inverse energy cascade was found in local box 



simulations as well ( Johansen et al 



2009). Here the turbulent advection term in the induction 
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equation drives large-scale radial magnetic field. 



The locality and anisotropy of the MRI turbulence is an important aspect for dust growth and 
therefor the planet formation. The eddies are stretched in the azi muthal direction d ue to the strong 



shear. They have a characteristic low tilt angle in the r -cp plane (iGuan et al. 



confirmed this ti 



Davis et al 



2OIOI : 



t angle for the velocity and the magnetic fields (IGuan et al. 



Guan & Gammie 



2011 



Sorathia et al 



correlation wavelengths is dependent on resolution (IGuan et al 



20091 : 



2009h. Several works 



Fromang 



2OIOI : 



2011b. The s ize of the corresponding 



20091) and conve rges by usin g 



a fixe d value of viscous and explicit dissipation in unstratified local simulations (|Fromang 



201 oh. Unstratified g lobal models interpret the magnetic tilt angle as con vergence parame ter 



( Sorathia et al. 



201 lb . They found convergence with tilt angles around 13°. 



Beckwith et al. 



(1201 Ih 



found tilt angles of 9° in global stratified simulations with spatial structures of the turbulent field 
in the order of H. 



Global disk simulat i ons (lArmitagd ll998l:lHawlev 



Fromang & Nelson 



201 ll: 



Sorathia et al 



2006, 



2009 



Dzvurkevich et al. 



2010 



2OOOI: lArlt & Rudiged 12001 



Flock et al 



2011 



Beckwith et al 



(2011 



( 2011 



2011b are used to study the MR I evolution on large s cales 



) found a stronger accretion stress compared to 



Beckwith et al, 



Fromang & NelsonI (120061) and 



Flock et al. 



) with a stronger initial toroidal field. Unstrati fied simulations sho w a similar correlation 



between accretion stress and the initial plasma beta (|Hawley et al 



1995b . Here a stronger seed 



field will drive to stronger accretion stress. The majority of stratified global disk simulations 
has been done for restricted (0 < njl) azimuthal domain sizes. At first glance, MRI turbulence 



behaves similar for both full 2n and smaller domain sizes ( Hawlev 



2000b . 



In our previous work we compared stratified simulati ons of nIA and 2;r in azimuth. There, we 



observe stronger az imuthal fields for the njA domain size (|Flock et al. 



global simulations ([Sorathia et al 



2011b . Recent unstratified 



2011b do not show large differences between domain size of 



n/A and In. This fact indicates a mean field dynamo mechanism. The stratification is crucial 
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for driving aCl dynamo i n disks and therefor for the creation of large scale magnetic fields 



(IKrause & Raedler 



1980|) . With this work we perform a detailed study of different azimuthal 



domain sizes. We investigate the turbulent and the mean field evolution for the velocity and 
magnetic fields. 

In stratified disk simulations, there is a periodic change of sign for the mean toroidal 
magnetic field. A similar periodicity of toroidal magnetic field, known as butterfly diagram, 
is observed in the sun. It could be explained by a MHD dynamo process. The MRI could b e 



self-sustaining by a analogous d ynamo process (|Hawley et al. 



Gressel 



2OIOI : 



Simon et al, 



19961 : 



Lesur & Qgilvie 



2008 



I; 



2011b . Strong shear in accretion disks will wind up any radial magnetic 



field generated by MRI and produce toroidal field. This field will act as seed for the MRI 



again. Solutions for aQ dyn amos in rotating systems were presented by 



1993h : 



Elstner et al 



ocal box simulations 



(120001) : 



Ruediger & Kichatinov 



19961). Calculations of the dynamo-a for MRI have been performed in 



Brand.enburg et al 



Ziepler & Rudiged (120001) 



dynamo-o;. iBrandenburg & Donned ( 



41995): 



Davis et al 



19971); 



Brandenburg & Donned (|1997|) : 



(201C): 



Gressel 



Rekowski et al 



(12(3101) showing a negatively 



Riidiger & PipinI (|2000l ) explained the negative sign 



as an effect of vertical b uoyancy. The first indications for a positive dyna mo-a were found in 



global disk simulations (|Arlt & Riidigei 



2001 



Arlt & Brandenburg 



200 lb . Dynamo solutions 



for positive or negative dynamo-a predict long-term global mean magnetic fields which become 
symmetric (quadrupole, dynamo-a"'"'''' < ) or asymmetric (dipole, dynarno-o;'""" ^^' > 0). E.g. 



dipole solutions support the creation of disk wind and lets (IRekowski et al 



of dynamo action in accretion disks was present e d by 



Brandenburg & von Rekowskil (l2007b : 



Blackman 



torn . 



200C). A review 



Brandenburg & SubramanianI (|2005h 



The connection between the dynamo processes and the large-scale magnetic field oscillations 



'Negative dynamo-ff means a negative correlation between the turbulent EMF (Electromotive 
force) and the mean toroidal field in the upper (northern) hemisphere. 
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was shown by 



Lesur & Ogilvid (l2008bh : 



Gressel 



are universal for stratified MRI simulations ( Stone et al 



(2010) 



Simon et al, 



1996; 



timescales of ten local orbits, presented rec ently in loc al (iGressel 



201 



Hawlev et al. 



201 



Flock et al 



Guan & Gammie 



2011 



Beckwith et al. 



201 ih a nd global dSorathia et al. 



( 2011 ). These oscillations 



Mi 



ler & Stone 



20 







200C) with 



Simon et al 



20101 : 



2011 



Dzvurkevich et al. 



20111) simulations. We use the second order Godunov 



code PLUTO which was successfully applied in recent global simulations (IFlock et al. 



2011 



Uribe et al 



2011 



Beckwith et al, 



2010L 



201 1|) . The paper is structured in the following way: 



First, we describe the disk model and the numerical parameter. For the results in section 3 we 
study the turbulent and the mean field evolution for all azimuthal domain. Section 4 and 5 present 
discussion and summary. 



2. Setup 



Our disk model is presented in detail in 
physical and numerical initial conditions. 



Flock et al. 



(|201l|) . We give here a summary of our 



Disk model 

The HD initial conditions of density, pressure and azimuthal velocity follow a hydrostatic 
equilibrium. We set 

/sin(0)-l\ 

with po = 1 -Oj H/R = Co = 0.07, R = rsin (6). The pressure follows locally an isothermal equation 
of state: P = c]p with Cs = Co/ V^. The azimuthal velocity is set to 

^ ^ r\ sin(e) 7 
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The initial velocities V,- and Vg are set to a white noise perturbation amplitude of Vl"f = 10 ^c^. 
We start the simulation with a pure toroidal magnetic seed field with constant plasma beta 
jS = 2P/B^ = 25. 

The radial domain extends from 1 to 10 AUg. The 9 domain covers + 4.3 disk scale heights, or 
6 = n/2 + 0.3. For the azimuthal domain we use four diff"erent models: (f)'^^'™' = ;r/4, n/2, n and 
2n. We use a uniform grid in spherical coordinates with an aspect ratio at 5 AU of 1 : 0.67 : 1 .74 
(Ar : r^e : rA0sin0). The resolution is fixed to A^,- : 384, Ng = 192 , A^^ = 768 • 0„,,„,/(2;r). We 
have around 23 grid cells per pressure scale height. 

Buffer zones extent from 1 to 2 AU as well as from 9 to 10 AU. In the buffer zones we use a 
linearly increasing resistivity to the boundary. This damps the magnetic field fluctuations and 
suppresses boundary interactions. In the buffer zones we use also a relaxation function which 
reestablishes gently the initial value of density over a time period of one local orbit. In the buffer 
zones we set: p"*^* = p - (p - p^"") • At/Toibits- Our outflow boundary condition projects the radial 
gradients in density, pressure and azimuthal velocity into the radial boundary and the vertical 
gradients in density and pressure at the 6 boundary. We ensure to have no inflow velocities. For 
an inward pointing velocity we mirror the values in the ghost cell to ensure no inward mass flux. 
The 9 boundary condition for the magnetic field are set to zero gradient, which approximates 
"force-free" - outflow conditions. The normal component of the magnetic field in the ghost cells 
is always set to have V ■ ^ = 0. 



^We set AU as unit length. As the simulations are scal e invariant, the rad ial extent could be 
also from 0.1 to 1 AU or from 10 to 100 AU, more details in Flock et all POIII ) 
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Numerical setup 



The detailed numerical configuration is pre sented in 



successfully used in recent global simulations by 



Flock et al. 



Beckwith et al. 



(1201 



(120101) and was also 



ll). For all runs we emplo y 



2005b, 



the second order scheme in PLUTO with the HLLD Riemann solver (iMivoshi & Kusand l! 
piece-wise linear reconstruction and 2'"' order Runge Kutta time integration. We treat the induction 
equation with the " Constrained Transport" (C T) method in combination with the upwind CT 



method described in 



Gardiner & Stonel (|2005|) . All models were performed on a Blue-gene/P 



cluster for in total over 3 million CPU hours. 



2.1. Measurement and integration 

For our analysis we use the central domair^ from 3 to 8 AU. Total volume integrations or a 
variable F, as used for the total stress are performed with 

ptotai ^ lFdV= Fr'^ sin edrdedc/). 

J J3 JOtegin Jo 

In global disk models, the gas dynamics are only self-similar along the azimuth. Therefor, mean 
values like v^, are always averaged over azimuth. This includes the calculation of the turbulent 
EMF' in Fig. 11. For further analysis we always use an 2D dataset of mean values, e.g. V^{r, 6) 
to construct the 3D turbulent dataset v^(r, 9, (p) = v^(r, 6, (f>) - v^(r, 6). For volume integration over 
mean values, as a'^^f", we use 



begin 

sin OdrdO. 



Some results are determined in the center of computational domain. Analysis done at 4.5 AU are 
the tilt angle calculations. Fig. 6, the mean field contour plots. Fig. 9, the parity. Fig. 10, the 



^The "central doma in" is here the dorn ain between 3 and 8 AU to avoid impact of the inner and 
outer buffer zones, (see lFlock et al.l (|201l|) ). 
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dynamo coefficients in Fig. 11. This results are averaged over azimuth and a small radial extent 
(+0.5 H = 0.16 AU). For the time evolution of the tilt angle, Fig. 6, top, we average vertically 
+0.5// at the midplane. Radial contour plots are averaged over azimuth and height, between 
- 1.5H. This applies for the mean toroidal field Fig. 3, the dynamo Fig. 1 1 and the mean fields 
in Fig. 12. The parity is averaged over the total disk height at 4.5 AU, Fig. 10. 



3. Results 



In this section we investigate the turbulent and mean field evolution for the azimuthal MRI for 
different azimuthal domain sizes. Table 1 summarises the results of accretion stress, contribution 
of mean magnetic field to the total stress, dynamo-a and RMS velocities for all models. Table 2 
summarises results of the two-point correlation function, including tilt angles, major and minor 
wavelength. For all models, the accretion disk becomes unstable to MRI on timescales of ten 
local orbits. All models develop an oscillating zero-net flux configuration after around 250 inner 
orbits. The time evolution of total magnetic energy. Fig. 1 left, is normalised over the total initial 
magnetic field energy B^. It shows the peak of magnetic energy shortly after the linear MRI phase 
around 100 inner orbits. Between 100 and 400 years, the total magnetic energy de creases due to 



loss of the net magnet ic flux and mass loss (see also Fig. 13 in 



Flock et al 



Beckwith et al 



(1201 ih and Fig. 3 in 



(|201 II) ). After 400 years, n/4 and n/2 models show strong fluctuations while n 
and In models do saturate. In the saturated state (> 800 inner orbits), the total magnetic energy 
evolution shows a relative constant level for the n and In model. 

All models have the same resolution per extent ((p^^^^"^/N^). The toroidal quality factor 
Qj) = Arrit/^(P shows the quality of res olved MRI (Q^ > 8). We follow the analysis done by 



(Noble etal 



20101 : 



Sorathia et al 



20111) and calculate the mean Qd, for the central domain (3 to 8 
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AU). The definition is similar to the toroidal quality factor 2^ by 



Hawlev et al 



( 2011 ) 



A(f) 



15 yj^* A(f) 



16 \B, 




15 ^/pQ.A(p 



Fig. 1, right, shows Q^p over time. For all models we have > 8. The n/4 and n/l show a higher 
due to stronger magnetic fields. 



3.1. Turbulent evolution - a value 

We start the comparison with the volume integrated turbulent stress scaled on the local 
pressure, e.g. the Shakura-Sunyaev ass ■ The ass value is determined from the turbulent Reynolds 



stress Tr = pv^v^ and Maxwell stress Tm = B^By^/Att. We split the total ass into a mean and 
turbulent component. For the Maxwell stress, we split the magnetic field components into the 
turbulent and mean component, e.g. = + B^. This leads to a second Maxwell stress 
component, e.g. the mean Maxwell stress T'^^'^ = B^ ■ BR/4;r. For the volume integrated turbulent 
ttg * value we integrate the mass weighted stresses over the central domain 



"ss - 



/pdV 

The same is done for the mean Maxwell stress 



mean 

"ss 



JpdV 



The volume integrated a'^g (Fig. 2 left - solid line) and the volume integrated or^l™ (Fig. 2 right 
- solid line) are plotted versus time. We are interested in the steady state and we use the time 
period between 800 and 1200 inner orbits for averaging. Fig. 2 (left) shows that the n/4 and n/2 
models present higher ass value than the n and In models. The mean magnetic fields provide a 
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Fig. 1. — Left: Total magnetic energy evolution ov 
time. All models show a well resolved MRI. 
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Fig. 2. — Left: Volume integrated 0:55 value for all models. Right: Volume integrated ass values 
using only the Maxwell component with the mean magnetic fields. Dotted lines show same results 
but for a 7r/4 average (0 - n/4) instead of whole domain size. 
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significant contribution to the total stress for the restricted azimuthal domains, see Fig. 2, right. 
The time averaged ratio between the turbulent Maxwell stresses and the mean Maxwell stresses is 
up to 33 % for the n/A model while it decreases in the full In model down to 8 %, see Table 1. In 
Table 1 we summarise the results of ass™, '^ss'' ^sT'- "^^^ standard deviation is determined by 
the temporal fluctuations. For model n/A we determine a^^^^^ = (1 1.8 + 2.3) ■ 10^^. For model n/l, 
ags^' reduces to (9.3 + 0.9) • 10 The stress of the two largest azimuthal domain sizes, n and In, 
matches within the standard deviation. For model n, the time averaged a^°^^ is (5.6 + 0.5) • lO""* 
and (5.4 + 0.4) • 10"^ for model In. 

To verify the results we made the same analysis in the same azimuthal extent for every 
model. Instead using the full azimuthal dataset for the analysis, we use here the azimuthal extent 
between - nIA in every model. The results are shown in Fig. 2, dotted lines. In Fig. 2, left, these 
ass values are only slightly lower than the total domain integration. This indicate that most of 
the turbulent stress is generated by the small scale turbulence (m < 8). In Fig. 2, right, these ass 
values represent the stress for one specific mode (m = 8). We see again that the smaller scales 
contribute more to the a^^^^ than the larger scales. We summarise that the turbulence is amplified 
in case for the nil and nIA model. These models present higher a^^^ and a™^^" values than the n 
and 2n runs. 



The n/A run presents another exceptional behaviour. Around 800 inner orbits, the a value 
increases quickly up to c = 0.013. The reason for this increase is connected to strong mean 
toroidal field oscillations. In Fig. 3 we plot contour lines of the resolved 1^"^^^ from the mean 
toroidal field 5^ with i^?/^*^ ^ 8. 



Accretion burst due to mean fields 
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The definition is equivalent to the definiton of the toroidal quality factor but calculated from 
the mean toroidal field instead from the total field (see Fig. 1, right). There is clear correlation 
between the rise of the ass value and resolved mean toroidal field. At the same time there is a 
superposition of strong mean field along radius, see Fig. 3 red solid line. The amplifications are 
present in the n/A model, Fig. 3 top, and the n/l model. Fig. 3 bottom. For the larger domains, n 
and 2n (Fig. 4), the mean field stays at lower values and Aj.^ is not resolved. 



Turbulent magnetic and velocity fields 

We investigate the spatial distribution of magnetic energy with Fourier analysis. The 
magnetic field amplitudes, ^jB{mY are plotted in Fourier space along azimuth at the midplane and 
for all models. Fig. 5, left. The plots show that the highest amplitudes of the magnetic fields are 
at the largest scales. The n/A and n/l model show systematically increased amplitudes compared 
to the n and In model. This is true for all modes and for all three magnetic field components. It 
is also visible in the time averaged total magnetic energy. Fig. 1 left dotted lines. Time averaged 
values, in units of the initial total magnetic energy, are B^/Bl = 0.54 + 0.12 for model n/A, 
0.48 + 0.09 for model n/2, 0.34 + 0.07 for model n and 0.35 + 0.07 for model 2n. Here, time 
average is done between 400 and 1200 inner orbits. We present the velocity field in Fourier space 
■\JV(my in Fig. 5, right. We observe increased turbulent velocities for the restricted domain 
models. The radial velocity (dashed line) dominates in the range between 2 < m < AO. The peak 
turbulent velocity is V,- at m = A for the n/2, n and 2n run. Coincidentally, this mode matches 
the domain size of n/2. The n/A does not include this mode. This lack of large scale turbulent 
radial fields becomes again vi sible in the velocity tilt ang l e. The peak at m = 4 (22H) is connected 



to spiral density waves. After 



Heinemann & Papaloizoul (|2009l) we should observe the peak at 



m = 14 (6H). This could be a resolution issue as the domain size of n/A (1 IH) should be large 
enough to include spiral density waves. 
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Fig. 3. — Contour lines of the resolved MRI from the mean toroidal field Aj^^ with the evolution 



of the a value for the models 7r/4 (top) and n/2 (bottom). The contour lines show A^^-^ 
strong mean toroidal field amplifies the turbulence. 
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Fig. 4. — Contour lines of the resolved MRI from the mean toroidal field Aj^^ with the evolution 
of the a value for the models n (top) and 2n (bottom). The contour lines show Aj^^ = 8. Here, the 
mean toroidal field is weaker and not resolved by the code. 
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Fig. 5. — Left: Magnetic field distribution in Fourier space over azimuthal wave number for all 
models and magnetic field components. Right: Same for the velocity field. Values are from the 
midplane and time averaged between 800 and 1200 inner orbits. 
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Parity 
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n/4 


11.8 


±2.3 


0.33 


8.9 = 


t 1.6 


-3.4 + 0.9 


3.3: 


tO.8 


-0.2 + 0.4 


0.125 + 


0.009 


n/2 


9.3 d 


bO.9 


0.19 


7.8 = 


tO.7 


-2.8 + 0.6 


3.1 : 


tO.7 


-0.2 + 0.5 


0.148 + 


0.006 


n 


5.6 d 


bO.5 


0.12 


5.0 = 


tO.4 


-2.4 + 0.3 


2.1 : 


tO.3 


-0.1+0.5 


0.112 + 


0.005 


2n 


5.4 d 


bO.4 


0.08 


5.0 = 


tO.3 


-2.3 + 0.2 


2.1 : 


tO.2 


0.2 + 0.4 


0.113 + 


0.005 



Table 1: Model overview. From left to right: Azimuthal domain; Volume integrated total stress; 
Relation between a^l™ to cfs*; a'^"* stress; Value of dynamo a^^ for southern hemisphere (lower 
disk); Value of dynamo a^f for northern hemisphere (upper disk). 
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Fig. 6. — Top left: Midplane magnetic tilt angle over time for all models. Bottom left: Time 
averaged magnetic tilt angle for all models. Top right: Midplane velocity tilt angle over time for 
all models. Bottom right: Time averaged velocity tilt angle for all models. 
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Two-point correlation function 



The two-point correlation function, specified for MRI by iGuan et al.l (120091) . allows to 
study the locality and anisotropy of the turbulence. We measure the tilt angle for the magnetic 
sin 2^5 = \B,B^\/B^ and the turbulent velocity field sin26'y = \V'^V'^\/V'^ at 4.5 AU. In Fig. 6, we 
plot the time evolution, top, and the vertical distribution, bottom, of the magnetic tilt angle 6b, left, 
and the velocity 6v, right. The time evolution of the magnetic tilt angle 9b is plotted in Fig. 6 top 
left. The n/4 and n/l model show higher tilt angles (Ob ~ 9°) with much higher time deviations as 
the n and 2n model (6b ~ 8°). The n/4 model shows sudden increase of the tilt angle at 80 local 
orbits. At this time, the turbulence gets amplified due to strong axisymmetric fields, see Fig. 3. 
The time averaged vertical profile of is plotted in Fig. 6, bottom left. The tilt angle present the 
highest values in the coronal region. Here, we see again higher 6b values for the n/4 and n/2. The 
n model shows smaller 6b at the midplane compared to 2n which is an artefact of the selected time 
average. Both models present equal values after 100 local orbits, see Fig. 6, top left. 

We do the same analysis for the velocity tilt angle dy. The time evolution for 6v does 
not show strong fluctuations. At the midplane, we measure a time averaged velocity tilt angle 
of 6v ~ 14° for all models except of n/4. The n/4 model shows a systematic lower tilt angle 

_ 22°. This becomes also visible in the vertical profile. Here all models, except n/4, show 
a peak of 6v at the midplane. The reason is unresolved density waves. The n/4 model does not 
resolve the density waves with m = 4. At m = 4, all models show the highest turbulent amplitude 
in the radial velocity. For model n/2 it matches the size of the domain and it is not captured by 
model n/4. The fast drop of magnetic and velocity tilt angles above 4 scale height could be due to 
boundary efi'ects. 

We calculate the two-point correlation functions in the r-^ plane: 6y =< 6Vi(x)6Vi(x + Ax) > 
and Eb =< 6Bi(x)6Bi(x + Ax) > with x = r, (/>. In Fig. 7 and Fig. 8 we present the two-point 
correlation function at 5 AU at 1 scale height with Ar = 2H = 0.1 AU and the total domain 
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rA0 = (f)^"'"""' /O.OIH. For the In model we have around 90H (In 10.01). The corresponding major 
and minor wavelength are calculated using the half width at half maximum (HWHM) in units of 
H (HI5AU = 0.35AU). It measures the distance be tween the center e =1.0 and e = 0.5 along the 



major Amuj and minor /l„,;„ axis, see footnote 7 in 



GuanetaL 



(|2009l) . We measure the two-point 



correlation function at different heights. The results between +2H are similar and we present the 
values at 1 scale height. For the velocity, the /Imaj of the nIA run is 1.1//. The n and In run present 
both a value of 1.9//. We find a similar increase for the /Imin, from 0.19// for nIA to 0.24// and 
0.23// for model n and In. The values of the 7r/2 model present the highest values, /Imaj = 2.0H 
and /Imin = 0.29H. This is again due to the peak of turbulent radial velocity at domain size, see 
Fig. 5, right. It is visible in the magnetic fields too. The /Imin value for the magnetic fields are 0.14 
H, except the njl model with 0.16 H. The A^a] increases with increasing the azimuthal domain, 
the nIA model with 1.1// to 1.4//, 1.6// and 1.7// for the full In. All resuhs of the tilt angels, 
major and minor wavelengths are summarised in Table 2. 

The models with nIA and njl show an amplified turbulence. The (p"'™' affects the large scale 
and small scale turbulent properties. Only an azimuthal domain of n does reproduce similar large 
scale and small scale turbulent properties as in the full In run. The strong mean field generated by 
the aO. dynamo are responsible for the MRI amplification. 



3.2. Mean field evolution 



A typical feature of MRI in stratified disks is an oscillating toroidal magnetic field, generated 
by oscillating radial magnetic field. This feature is well known as 'butterfly' pattern, which wings 
appear due to the buoyant movement of the toroidal field from the midplane to upper layers. The 
timescale of these oscillation is around ten local orbits. Recent work in local box simulations 
showed the context between this osci 



Simon et al. 



2011 



Hawlev et al, 



2011 



l ating magnetic fie 



Guan & Gammie 



d and a dynamo process (|Gressel 



20101 : 



201 1|) . In this section we investigate the 
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n/A 12.0 1.1 H 0.19H 9.1 1.1 H 0.14H 

n/l 14.1 2.0 H 0.29 H 8.9 1.4 H 0.16 H 

n 14.1 1.9 H 0.24 H 7.7 1.6 H 0.14 H 

In 14.2 1.9 H 0.23 H 8.2 1.7 H 0.14 H 



Table 2: Two-point correlation values for all runs. From left to right: Azimuthal domain, correla- 
tion angle for the velocity, wavelength of the major axis, wavelength of the minor axis, correlation 
angle for the magnetic field, wavelength of the major axis, wavelength of the minor axis. 
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Fig. 7. — Contour plot of the two-point velocity correlation function at 1 scale height at 5 AU. The 
red line shows zero contour. 



-21 - 




Azimuthal [H] Azimuthal [H] 

Fig. 8. — Contour plot of the two-point magnetic field correlation function at 1 scale height at 5 
AU. The red line shows zero contour. 
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evolution of this axisymmetric magnetic fields and the connection to the dynamo process. 

The parity and butterfly pattern 

In Fig. 9, top, we present the time evolution of axisymmetric radial and toroidal magnetic 
field over height. The values are normalised over the initial toroidal field. The generated toroidal 
magnetic field, Fig. 9 (second from top) is around one order of magnitude higher than the radial 
magnetic field. We observe a change of sign every 5 local orbits. The butterfly wings are mostly 
antisymmetric with respect to the midplane. To quantify the symmetry we determine the parity 
of the mean magnetic field. We calculate the symmetric (S) and asymmetric (AS) magnetic field 
component: 5^^^ = 0.5(8^^^ + 5^^^) and Bf^^ = 0.5(8^^^ - 5^^^) with the values of the northern 
(NH) and southern (SH) hemisphere (Snfl The parity 

is determined with total dipole and quadrupole energy components E° = (B^^)^ + (B^)^ + (B^^)^ 
and = (Bf )^ + (B^^)^ + (B^)^. The toroidal field is much larger then the radial and theta 
magnetic field. It is possible to define a symmetric (Quadrupole) or antisymmetric (Dipole) 
configuration as the total parity is set by the toroidal field. Then, a parity of -1 defines a 
pure symmetric configuration (Quadrupole) while a parity of +1 defines a pure antisymmetric 
configuration (Dipole). The time evolution of the total parity is plotted in Fig. 10, top, for all 
models. The total parity starts with -1 as the initial field is symmetric. The parity of only B, 
and Bq is plotted in Fig. 10, bottom, and present a similar time evolution. Both parities change 
sign several times during the simulation for all models. The time averaged values (400 - 1200 

"^The northern hemisphere is placed on the upper disk if the azimuthal velocity is positive. Then 
if one looks at the north pole , the disk is rotation counter-clockwise in the northern hemisphere, 
e.g. mathematically positive. 



inner orbits) show strong deviations around zero parity, see Table 1 . Only the 2n model is mostly 
antisymmetric for the simulation time. The contour plot of total parity over height, Fig. 9 third 
plot from top, shows the correlation between the parity and the 'butterfly' pattern. The symmetry 
of the mean toroidal field in respect to the midplane sets the total parity. Even the total parity is 
mostly antisymmetric (yellow, +1) there is a change of the parity to symmetric for two butterfly 
cycles between 80 and 100 local orbits (also visible in Fig. 10, solid line). 



orQ Dynamo 



In mean field theory, there is a mech anism to generate 



turbulent field. In case of an aO. dynamo (IKrause & Raedlei 



large- scale magnetic fields by a 



19801) there should be a correlation 



between the turbulent toroidal electromotive force (EMF'^) component and the mean toroidal 
magnetic field. 



EMF^ = a^^B^ + higher derivatives of B 

with EMF'^ = v'.5g - v'gB'^. The sign of a^^ has to change for the southern and northern 
hemisphere. The correlation is plotted in Fig. 11, left, for the northern hemisphere (top) and the 
southern hemisphere (bottom). We get a positive sign for the a^^ in the northern hemisphere (a^f) 
of the disk (Fig. 1 1 top) and a negativ e sign in the southern hemisph ere (ff?^). This result was 



predicted fo r stratified accretion di sks (|Ruediger & Kichatinov 



simulations (|Arlt & Rudiger 



19931) and also indicated in global 



20011) . Each dot in Fig. 1 1 left, represent a result from a single time 



snapshot. The boxes show the limits of the values for each model. The n/A- and n/2 model show 
higher amplitudes in the mean field as well as in the EMF'^ fluctuations. All values of a^^ 
are determined using a robust regression method and summarized in Table 1 . A time evolution 
of the mean field and the turbulent EMF'^ is presented in Fig. 11, right, for model In, top, and 
model 7r/4, bottom. In Fig. 1 1 right, we divide the turbulent EMF^ with the measured (see 
also Table 1). The n/A run shows higher fluctuations compared to the In run. A time evolution 
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Fig. 9. — Top to bottom: 1. Mean radial magnetic field over height and time. 2. Mean toroidal 
magnetic field over height and time. 3. Contour plot of the parity over height and time. 4. Contour 
plot of B^a^^ I EMF'. over height and time. All plots are made for model In at 4.5 AU. 



-25- 



Dipole dominated (Parity > 0, antisymmetric) 

7C/4 




20 40 60 80 100 120 

Local orbits 




Fig. 10. — Parity of mean toroidal magnetic field (top) and of mean poloidal field (bottom). 
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Fig. 11. — Top left: Correlation between the mean toroidal magnetic field and the turbulent EMF 
component EMF'^ for the northern (upper) hemisphere of the disk and for all models. Rectangles 
show the limits of the data values. Bottom left: Correlation between the mean toroidal magnetic 
field and the turbulent EMF component EMF'^ for the southern hemisphere of the disk and for 
all models. Top right: Time evolution of mean toroidal field (solid line), over-plotted with the 
turbulent EMF (red dotted line) divided by a™ for model In. Bottom right: Time evolution of 
mean toroidal field (solid line), over-plotted with the turbulent EMF'^ (red dotted line) divided by 
a^Jf for model njA. 
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of ■ cf^^/EMF^ over height is presented in Fig. 9, bottom. We see that the sign of is well 
defined for the two hemispheres, reaching up to 3 scale heights of the disk. 

3.3. Mean fields over radius 

In this section we study the development of the mean magnetic fields along radius. We show 
results from our full In model as it represents the most realistic physical domain size. A contour 
plot of mean toroidal field, normalised over the square root of the pressure, is presented in Fig. 
12, top right, over radius and time. All results in Fig. 12 are averaged along azimuth and along 
d between the midplane and two disk scale heights in the northern hemisphere. Fig. 12, top 
right, shows the irregular change of sign for the mean toroidal magnetic field along radius. The 
timescale of the "butterfly" oscillations at a given radius can change because of radial interactions. 
The timescale of reversals of the toroidal magnetic field does vary from the ten local orbital line 
(see Fig. 12, top right, horizontal homogeneous B^. The mean field configuration along radius 
can strongly afi'ect the accretion stress, see Fig. 3. The distribution of mean over radius is more 
irregular compared to the toroidal field, see Fig. 12 bottom left, although we observe a preferred 
sign of mean Bq for a specific radial location, e.g. positive over time between 4 and 5 AU. A 
time evolution over radius of ■ a^^/EMF^ , Fig. 12 top left, shows again the positive sign 
of a^^ in the northern hemisphere (see also Fig. 9, bottom). By definition, the a^^ presents the 
same distribution along radius as the mean toroidal magnetic field. In contrast we do not find a 
correlation between the turbulent velocity of the gas and the distribution of mean magnetic fields. 
Fig. 12, bottom right, presents Vnns over radius and time for the northern hemisphere. The RMS 
velocity is about 0.1 c.,, nearly constant over radius and time. A time average of Vni,s is given in 
Table 1 for all models. We again emphasize the lower turbulent velocity in the ;r/4, compared to 
7r/2, due to the lack of the radial velocity peak (see Fig. 5, right). 

In our previous work we have shown the 1 fr profile for the turbulent magnetic fields 
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Fig. 12. — Top left: Contour plot of ■ a^^ /EMF'^ over radius and time (see also Fig 9, bottom). 
Top right: Mean toroidal magnetic field over radius and time. Bottom left: Mean 6 magnetic field 
over radius and time. Bottom right: Turbulent RMS velocity over radius and time. All plots are 

marlp. fnr mnrlp.l Ott in thp. nnrthp.m hp.misnhp.rp. 
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Fig. 13. — Radial distribution of the peaks of mean toroidal magnetic field. Values from the 
northern hemisphere are in red and from the southern hemisphere in blue, (see also Fig. 12, top 
right). 



-30- 



(IFlock et al, 



201 II) . Because of the time oscillations, it is difficult to estimate a radial profile for 
the mean magnetic field. To determine a time averaged radial profile of the mean toroidal field we 
measure the amplitude values of the oscillations. We use five different radial locations to measure 
the peak values of the mean toroidal field. The results are plotted in Fig. 13 for the southern (blue) 
and northern hemisphere (red). The amplitudes of mean toroidal field decreases with radius. 
The relative low number of values and their high standard deviation makes it difficult to fit. A 
1/r profile would apply (Fig. 13, green solid line). The values in both hemispheres look quite 
symmetric (Fig. 9, blue and red) and we do not see a preferred hemisphere for the mean field 
generation. 



4. Discussion 



After the saturation of MRI, the initial magnetic field configuration is lost. Each model 
develops oscillating mean magnetic fields which appears to be strongest in the n/A and n/l run. 
The strength of turbulence follows this trend. The mean fields are generated by a dynamo process 
which relies on the symmetry and on the strength of the turbulent field. We measure higher 
dynamo coefficient for the n/l and njA model as well as higher Maxwell st resses. This agre es 



wiffithe correlation between Maxwell stress and dynamo coefficient, found by 



Rekowski et al 



(|2000|) . The effect of increased magnetic energy at domain size seems to be indepen dent of 



resolution in stratified simulations (compare Fig. 12, bottom left, mode 
(1201 1|) ) but not present in unstratified simulations (compare Fig. 9b in 
they do not develop a dynamo. 



FO and PO in 



Sorathia et al, 



Flock et al 



(1201 ih ) as 
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Energy pile up and magnetic dynamo 



Which physical process is sensitive to the domain size and lead to the increased mean toroidal 
fields in n/A and 7r/2 models ? The first mechanism leads to the dynamo process as it generates 
axisymmetric magnetic fields out of the turbulence. Ano ther way to transport m agnetic energy at 



domain size could be due to an inverse energy cascade. 



Johansen et al 



(120091) showed in local 



box simulations that the Keplerian advection term in the induction equatio ns drives an inverse 



Riidiger et al, 



(l2007h 



energy cascade. This will lead to a transport of energy to larger scales. Also 
found in Taylor-Couette experiments that MRI, launched from a toroidal field, will have most of 
magnetic energy at the m = 1 and m = mode. 

Another open question is the sign of in global simulations. We find a positive o-^^, 
independent o f the azimuthal domain size. This positive has been indica ted for global 



simu 



1995 



2010 



ations by 



Brandenburg & Donnei 



Gresse] 



1997 



Arlt & Rudigei (1200 lb. Local simulations show a negative (Brandenburg et al 



Rudiger & Pipin 



200ol : 



Ziegler & Rudigei 



200ol : 



Davis et al. 



20101) . The reason of stronger mean fields in reduced azimuthal models as well as 
the positive sign of in global simulations has to be investigated in future work. One possibility 
would be to implement the 'Test field' method and to measure other components of the dynamo 



and diffusivity tensor, as it was done in 



Gressell ( 20101) . 



4.1. Time variability of accretion stress. 

Oscillating mean field are organized in elongated radial patches, normally following the 
time-line of ten local orbits. It can occur that for a given time, mean toroidal field of one sign 
covers the whole radial extent (3-8 AU). In such a case, temporal linear MRI will lead to a 
peak in accretion stress. Fig. 3. The effect of mean toroidal field, stretching over the whole 



radius, is independent on the azimuthal domain size, compare Fig. 12 top right. The amplification 
of accretion stress due to linear MRI, is visible only in the n/A model, as it present strongest 
amplitudes in the mean toroidal magnetic field. 



Correlation functions 



Beckwith et al, 



( 2011 ). We 



We confirm the results of recent stratified global simulations by 
find similar correlation angles (around 9°) and wavelengths (around H) for the magnetic field. 
A larger correlation length is expected because of the relative low resolution per scale h eight 



compared to local simulations (iGuan et al 



unstratified global simulations 



2009; 



Sorathia et al. 



Hawlev et al. 



2011 



Sorathia et al. 



2011). Recent 



(|2011l) suggest a magnetic tilt angle of around 13" 
for converged MRI turbulence. It remains still unclear how this could be applied for stratified 
disks with a minimum of 6b at the mi dplane. We found a magnetic tilt angle of around 13° 



above 2 scale heights. As discussed in 



Flock et al. 



(|201lh we believe to find convergence with 



resoluti ons around 32/64 gri d cells per pressure scale height. Here, a Fargo MHD approach as 



used in 



Sorathia et al. 



mm would be helpful. 



5. Summary 

We have studied the impact of different azimuthal extents in 3D global stratified MHD 
simulations of accretion disks onto the saturation level of MRI with an initial toroidal magnetic 
field. 

• Turbulence in restricted domain sizes like n/l and nIA is amplified due to strong toroidal 
mean field oscillations. For these runs, the Ac-it of the mean field is resolved leading to a 
temporal magnification of the ass value and increased total magnetic energy. In addition, 
radial superpositions of such strong mean fields can drive to a strong episodic increase of 
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accretion. The time averaged total ass is 1.2 + 0.2 ■ 10 ^ for model n/4, 9.3 ± 0.9 ■ 10 ^ for 
model n/2 and converge to 5.5 + 0.5 • 10"^ for both models n and 2n. 

• We find a positive dynamo a^^p for all models, a positive correlation between the turbulent 
EMF'^ and the mean toroidal magnetic field in the upper (northern) hemisphere. For the In 
model we found or^"'*^ = 2.1 + 0.2 • 10"^. The n/l and n/A present higher a^^ values but 
with stronger fluctuations in EMF'^ and mean B^. 

• The 7r/4 and 7r/2 models show higher tilt angles and smaller correlation wavelengths in 
the two-point correlation of velocity and magnetic field compared to the n and In models. 
We find e^^^ = 14° for models > n/2 and 0™' = 12° for model n/4. The n/4 model does 
not resolve the peak radial velocity at m = 4. The tilt angles for the magnetic fields are 
smaller. At the midplane we observe time averaged magnetic tilt angles between = 8 - 9° 
increasing up to 9b = 12 - 13° in the corona. For the full 2n model we found A'^^ = 1.9H 
andi-J = 1.7H. 

• The parity of the mean magnetic fields is a mixture of dipole and quadrupole for all models. 
The total parity is set by the oscillating toroidal field. The timescale of symmetry change 
between dipole and quadrupole is around 40 local orbits. The time evolution of the parity is 
distinct in each model. The 2n model remains longer in a dipole (antisymmetric) dominated 
configuration for the simulation time. 

We conclude: In global MRI simulations of accretion disks an azimuthal domain of at least 
n (180°) is needed to present the most realistic turbulent and mean field evolution as the full 
2n model. Here, the aO. dynamo plays a key role in determining the saturation level of MRI. 
Restricted domains of n/4 and n/2 amplify the MRI turbulence due to a stronger axisymmetric 
magnetic fields. 
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